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Large scale galactic turbulence: can self-gravity drive the 
observed HI velocity dispersions? 
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ABSTRACT 

Observations of turbulent velocity dispersions in the HI component of galactic discs 
show a characteristic floor in galaxies with low star formation rates and within indi- 
vidual galaxies the dispersion profiles decline with radius. We carry out several high 
resolution adaptive mesh simulations of gaseous discs embedded within dark matter 
haloes to explore the roles of cooling, star-formation, feedback, shearing motions and 
baryon fraction in driving turbulent motions. In all simulations the disc slowly cools 
until gravitational and thermal instabilities give rise to a multiphase medium in which 
a large population of dense self-gravitating cold clouds are embedded within a warm 
gaseous phase that forms through shock heating. The diffuse gas is highly turbulent 
and is an outcome of large scale driving of global non-axisymmetric modes as well as 
cloud-cloud tidal interactions and merging. At low star-formation rates these processes 
alone can explain the observed HI velocity dispersion profiles and the characteristic 
value of ~ lOkms"^ observed within a wide range of disc galaxies. Supernovae feed- 
back creates a significant hot gaseous phase and is an important driver of turbulence 
in galaxies with a star- formation rate per unit area > 10~^ yr~^ kpc~^. 
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1 INTRODUCTION 

The interstellar medium (ISM) is do minated by ir- 
regular/turbulent js^as motions (e.g. iLarso 3 Il98ll : 
lElm egreen Scald 120041 ). HI emission lines in most 
spiral galaxies have characteristic velocity dispersions of 
a ^ lOkm/s on a scale of a few hundred parsecs, exceeding 
the values expected fr om purely thermal effects. The data 
in Fig.[Tl assembled bv lPib et al.l (^006), also shows a tran- 
sition to much larger values in active/starbursting galaxies. 
Recen t high resolution observations by Petric & Rupen 
(|2007|1 of the nearly fa ce on disc galaxy NGC 1058 (see 
also Ibickev et al.1 Il990h provides us data on the radial 
behavior of the vertical velocity dispersion. They find that 
the dispersion declines with radius from ^12 — 15 kms~^ 
in the inner parts to ~ 4 — 6 kms~^ in the outer and is 
uncorrelated with active regions such as star formation 
sites and spiral arms. This is attributed to small scale 
(< 0.7 kpc) bulk motions. Petric & Rupen state that any 
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model attempting to explain turbulence in the ISM must 
also explain the radial decline that also has been detected in 



previous studies of e.g. NGC 6946 (Boulanger & Viallefond 
1^), NGC 628 ( Kamphuis Sancis i 1993; van der Hulst 
IS), NGC 2915 iMeurer et aLlll996l ). 

The main source (s ) of e nergy driving the ISM dynam- 
ics is still not clear ([Burke rt 2006), even though there 
are se veral candidates capable o f driving the ISM turbu- 
lence (|Mac Low KlessenI 120041 ) . A commonly discussed 
source is of stellar origin i.e. large -scale expanding outflows 
from high-pressure HII regions teessel-Dev net BurkertI 
r2003), stellar winds or supernovae. Whilst supernovae ex- 
plosio ns might dominate the energy input i nto the ISM (e.g. 
Mac Low &: KlessenI [iooi : IPib et al.|[2006h . the mechanism 
is unable to explain the broad HI lines in galaxies with 
a low star formation rate (SFR) and in regions of mod- 
erate stellar activity as in the outer parts of disc galax- 
ies. Many numerical studies have been carried out to un- 
derstand the infl u ence of supernovae in galac t ic dis c s (e.g 



Kim et al. I2OOII : Ide Avillez Breitschwerdtl l2004l. l2005l: 



Slvz et al.f I2OO5I : iMac Low et al.1 l2005l : I Joung k Mac Lowl 
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I2OO6I I. [Pib et all (|2006h reproduced the starbursting tran- 
sition seen in Fig. [1] but was unable to produce velocity dis- 
persions larger than ~ 3kms~^ for low values of SFR/Area. 
This strongly suggests that something else is contribut- 
ing to the energy budget. In addition, large scale holes, 
usually attributed to correlated supernovae explosions (e.g. 
IPuche et aLlllOQ^ ). are in some cases surp risingly uncorre- 
lated to stellar activity ([Rhode et al.lll999h . 

Another source of turbulence is galactic rotation. This is 
a huge reservoir of energy (Fleck 1981) and any mechanism 
able to generate random motions from ordered circular mo- 
tion could s ustain turbulence f or m any orbital times. Numer - 
ical work of IWada et all (|2Q02l ) and lWada k NormanI (|2007l ) 
has shown that realistic global models of galactic discs 
form a very complicated turbulent velocity field associated 
with a multiphase ISM. The only active source for this is 
shear coupled to gravitational and thermal instability. Local 
isothe rmal simulations of the I SM done by Kini Ostrikerl 
(I2OO7I) fals o previous work e.g. iKim Ostrikerl feOQl) and 
Kim et al. 1 f2003)) support this notion. They demonstrated 
that gas in a marginally stable galactic discs obtains, un- 
der certain conditions, velocity dispersions as large as the 
sound speed (here Cg — 7kms~^ ) due the swing- amplifier 
iGoldreich k Lvnden-Belll Il965bl : ljulian k Toomrd Il966l : 
iToom re 1981 ; Fuchs 2001 j. The swing- amplifier is when a 
leading wave is amplified into a trailing wave. The underly- 
ing mechanism is shear and self-gravity. 

iFukunaga k Tosal (Il98£ 



showed that rotational energy 
randomizes the motions of the cold cloud component of 
a galactic disc via gravitational scattering from their ran- 
dom epicyclic motions. This was was later quantified by 
iGammie et al.l (jT991|) who showed that the cloud velocity 
dispersion could reach ~ 5 — 6 km s~^ in this way, in agree- 
ment with observations (jStark k Brandlll989l ). We will dis- 
cuss this mechanism and its impact on the ISM in more 
detail in Sect|3l 

The Magneto-Rotational-Instabilit y (MRI) 

(|Balbus Hawlevlll99ll : ISellwood k Balbuslll99 9') coupled 
with galactic shear i s also a po ssible driver of turbulence . 
IPiontek k Ostrikerl (|2004l ) and IPiontek k Ostrikerl (|2005l ) 
obtained reasonable values of ~ 8kms~^ under favorable 
conditions. This mechanism becomes significant at low 
densities and might be important in the more diffuse outer 
part of galaxies. 

In this paper we carry out high-resolution 3-dimensional 
Adaptive Mesh Refinement (AMR) simulations to form a re- 
alistic multiphase ISM in which we can disentangle the con- 
tributing effects of self-gravity and supernovae driven tur- 
bulence. The simulations incorporate realistic prescriptions 
for cooling, star formation and supernovae feedback. Simi- 
lar numerical sin i ulatio n s have been carried out before (e.g . 
Gerritsen k Icke' "1997'; 'Wada et al." "2002*; iBottemal l2003l : 
I^sker k Brvan 2006; Wada k Norman 200l) but without 
addressing directly the issues discussed in this paper. 

The paper is organized as follows. In Sect. [2] we describe 
the numerical method used for this work and the setup of the 
galactic discs. In Sect. [3] we present the results from the nu- 
merical simulations, where the results treating the turbulent 
ISM are given in Sect 13.31 Sect. |4] summarizes and discusses 
our conclusions. 
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Figure 1. Characteristic HI gas velocity dispersion of a sample 
of galaxies as a function of the derived star formation rate in 
units M0yr-ikpc~2 as plotted by Dib et al. (2006). The figure 
is reproduced here by courtesy of Sami Dib, Eric Bell and Andreas 
Burkert, and by permission of the AAS. 



2 NUMERICAL MODELLING 

2.1 The code and subgrid modelling 

We use the adaptive mesh refineme nt (AMR) hydrodynam- 
ics code RAMSES (jTevssied [20o3 l . The code uses a sec- 
ond order Godunov scheme to solve the Euler equations. 
The equation of state of the gas is that of a perfect mono- 
atomic gas with an adiabatic index 7 = 5/3. Self- gravity of 
the gas is calculated b y solving the Poisson equation using 
the multigrid method terandtl llQTTl l on the coarse grid and 
by the conjugate gradient method on finer ones. The col- 
lisionless star particles are evolved using the particle-mesh 
technique. The dark matter is treated as a smooth back- 
ground density field that is added as a static source term 
in the Poisson solv er. The cod e adopts the cooling function 
of Sutherland k D opital (| 19931 ) for cooling at temperatures 
10^ - 10^-^ K. We extend coohng down to 3 00 K using the 
parametrization of [Rosen k Bregmanl (|l995l l. The effect of 
metallicity is approximated by using a linear scaling of the 
functions. 

The star fo rmati on recipe is described in 
[Dubois k Tevssierl (|2008l l but we summarize the main 
points here for completeness. In a cell, gas is converted to a 
star particle using a Schmidt law 



-— ifp > po 

6* 



= otherwise, 



(1) 



where t* is the star formation time scale and po is an ar- 
bitrary threshold that should be chosen to carefully make 
physical sense when related to the resolution and cooling 
floor. The star formation timescale is related to the local 
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free-fall time, 



-1/2 



(2) 



The parameters po and to are in reality scale dependent and 
not very well understood theoretically. A common way to get 
around this is to calibrat e them t o star formation rates in 
local galaxies i.e. to the KennicuttI 0-998) law and make sure 
that the values are cor npatible with modern e stimates of star 
formation efficiencies (jKrumholz Tanl[2007. ^ of 1 - 2 % 
per free-fall time in giant molecular clouds (GMCs). For 
example, if the star formation threshold po = 100 cm~^, the 
free fall time is 5 Myr meaning we can use to = 250 Myr 
to get 2 % efficiency per free fall time. As soon as a cell 
is eligible for star formation, particles are spawned using a 
Poisson process where the stellar mass i s connected to the 
chose n threshold and code resolution fsee lPubois Tevssierl 
l2008 l). 

The implementation of supernovae feedback is also de- 
scribed in the above reference (see their Appendix A). In 
the simulations that include feedback we assume that 50 % 
of the total supernovae energy, Esn = 10^^ ergs, goes into 
thermal energy where ?7sn = 10 % of each solar mass of stars 
that is formed is recycled as supernovae ejecta. The energy 
and gas release is also delayed by 10 Myr from the time of 
explosion by creating debris particles on at the time of explo- 
sion. By delaying the energy and mass release we allow for 
it to take place outside of dense environments, hence pre- 
venting it fro m radiating away too quic kly. We follow the 
prescription of Dubois & Teyssier (2008) and apply a large 
mass loading factor for the debris particles, i.e. r]w = 1-0. 
We use thermal rather than kinetic energy releases since we 
resolve the clumpy ISM and follow shocks self-consistently. 
In addition, a model that allows for debris particles to trans- 
fer kinetic energy is no longer valid as a Sedov explosion 
assumes a homogeneous medium to propagate into. While 
any treatment is inherently sub-grid, the supernovae impact 
should converge with enough resolution ([Ceverino Klvp inI 
I2OO7I I. 

In order to model subgrid gaseous equation of states and 
to avoid artificial gas fragmentation, the gas is given a poly- 
tropic equation of state as it crosses po. The temperature is 
set to 



where ro is the scale radius and h{r) the scale height. The 
sech^-term owes to the isothermality of the gas such that 

where Cs is the local sound speed and S(r) the local total 
surface density (S = Sgas + Sdm), naturally leading to a 
flaring disc. Experiments using a radial 1/r gas distribution 
were also performed (not discussed further) without any sig- 
nificant differences. 

We choose to model an M33 type galactic disc as it is a 
nearby well observed gas rich system. All global characteris- 
tics of the ini tial disc are in a greement with the observations 
presented bv lCorbellil (|2003[ ). M33 has a total gas mass (HI 
+ HII + He) of ~ 3.2 X lO^M© and an estimated stellar mass 
of 3 — 6 X 10^ Mq. We choose the initial gas mass to be in 
the high end of the total baryonic mass, i.e. ~ 9.2 x lO^M© 
as a lot of the outer material will not be a part of what we 
would associate with a galaxy and is only an artifact of the 
way we model isolated discs. The gas is assumed to have a 
mean metallicity of 0.3^0. 

We initialize the disc in a stable configuration where 
most of the disc has a Toomre parameter Q ^ 2 — 3. These 
large values of Q are desirable as we want the cooling to 
initiate instabilities and not our choice of initial conditions. 
A fairly large scale radius of r = 4.0 kpc is used. As this only 
reflects the very early setup of a forming disc galaxy this 
will not be an issue in our modelling. The dark matter halo 
has a concentration of c = 8.0, scale radius Rs = 35.0 kpc 
and a total mass of lO^^M©. These model parameters can 
be perceived as odd but is necessary for a best fit NFW- 
halo which, in accordance with observations (Co rbe lli 2003 ), 
reproduce a dark halo mass that within 17 kpc is ~ 5 x 
1O^°M0. This mass sets a lower limit on the actual dark 
matter halo mass. The HI velocity profile is still (slowly) 
rising at this radius. Scenario with different mass profiles, 
gas masses, shear and cooling floors are also explored. 

Numerically this setup is initialized at a resolution of 
100 pc using a nested hierarchy of grids situated in a simula- 
tions cube of size Lbox = 200 kpc. We achieve higher resolu- 
tion by refining cells both on based on a density and Jeans 
mass criterion (see Sect. 12.41 for details). The maximum al- 
lowed resolution is indicated in Table 1. 



To ( ^ 

Po 



(3) 



To is set to be the cooling floor of our simulations for con- 
sistency and 70 = 2.0. 



2.2 Initial conditions 

Our initial condition (IC) is an axisymmetric g alactic gas 
disc in equilibrium with an NFW ([Navarro et al. l 1997) dark 
matter halo. All relevant IC characteristics are presented in 
Fig. (|2]). The disc is initially isothermal at T = 10^ K having 
an exponential density profile, in cylindrical coordinates r 
and z, 



p{r,z) — sech^{z/h{r))poe 



(4) 



2.3 Simulation suite 

The performed runs are listed in Table 1. RUNO serves as our 
base run where we only consider the self- gravitating cooling 
gas and dark matter. RUNl introduces star formation as 
does RUN2 but at a higher resolution. RUN3 is identical to 
RUNl but implements the feedback prescription described 
m Sect.ET] These are our four fiducial simulations to un- 
derstand the importance of these physical mechanisms. To 
explore how choices of the gas cooling changes the outcome, 
RUN4, 5 and 7 are identical to RUNO except for a trunca- 
tion in the cooling function at the indicated thresholds. This 
will determine the ab ility of a disc t o develop a gravitotur- 
bulent state fsee e.g. Gammi 3l200ll). RUNG adopts 1/3 of 
the gas mass, making it a much more stable system. Finally, 
RUN8 adopts a very concentrated halo (c = 40, Rs = 7 kpc, 
M = 3 X lO^^Mo) peaking at 120 kms"^ to assess the influ- 
ence of a different shear, dQ/dR. 
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Figure 2. Characteristics of the initial gas distribution in all simulations, apart from RUNG and RUNS. The panels show, from left to 
right, the radial dependence of the surface density, the gaseous Q- value, the circular velocity and the scale height. 



Table 1. Performed simulations. 



Simulation 


Min. Ax 


Po [cm 3] 


Cooling floor, Tq 


Modelling comments 


RUNO 


24 pc 




300 K 


Hydrodynamics with self-gravity 


RUNl 


24 pc 


10 


300 K 


Like RUNO + star formation 


RUN2 


6 pc 


100 


300 K 


Like RUNl 


RUN3 


24 pc 


10 


300 K 


Like RUNl + supernova feedback 


RUN4 


24 pc 




1000 K 


Like RUNO 


RUN5 


24 pc 




10 000 K 


Like RUNO 


RUN6 


24 pc 


10 


300 K 


Like RUNl but with 1/3 of gas mass 


RUNT 


24 pc 




5 000 K 


Like RUNO 


RUNS 


24 pc 


10 


300 K 


Like RUNl but higher halo concentration 



2.4 Numerical considerations 

In these types of experiments it is important to consistently 
resol ve the Jeans sc a le ass ociated with the chosen cooling 
floor. iTruelove et aD (|l997l l demonstrated, using isothermal 
simulations, that at least 4 resolution elements are neces- 
sary to avoid artificial gas fragmentation. Our strategy is to 
choose realistic star formation density thresholds that to- 
gether with the cooUng temperature floor gives us a Jeans 
scale that can be resolved according to the Truelove criteria. 
The Jeans length given by 



(6) 




In our simulations, where the temperature floor is set at 
T — 300 K, this can be rewritten as 

(7) 

where n is expressed in cm~^. In order to satisfy the 
Truelove criterion we adaptively refine on a Jeans mass 
down to a resolution of Ax = 24(6) pc in RUNl (RUN2) 
where we have set the star formation density threshold to 
pQ — 10(100) cm~^. In addition to this precaution, the back- 
ground ISM polytropic EOS (see Sect. 12.1)) is activated at 
the same threshold, ensuring us that the Jeans scale never 
falls below the minimum value of - 100(25) pc set by Eq.[Zl 
(at n = 10(100) cm~^). Additional simulations have been 
conducted adopting 16 cells per Jeans length to asses the 
fidelity of the gravitational fragmentation without any sig- 
nificant difference in outcome. No physical perturbations for 
gravitational instability are seeded in the initially smooth 
disc meaning the actual perturbations existing arise from 
the AMR grid. As we are only interested in the long-time 



(t > 1 Gyr) dynamical evolution of the system, the actual 
morphology of the early unstable disc is of little importance. 



3 RESULTS 

As we will describe in Sect. 13.31 the response to gravita- 
tional instabilities, both in the form of bound structures 
and local non-axissymmetric instabilities, is an important 
source of turbulence. Therefore, before addressing the is- 
sue of turbulence we characterize the global gas evolution, 
phase-structure (density and temperature) and stability of 
the simulated galactic discs in Sect. l3.ll and Sect. 13.21 



3.1 General evolution and morphology 

3.1.1 Gas evolution 

Fig. [3] shows the total gas surface density maps of RUNl, 
RUN2 and RUNS at t = 0.5, 1.0, 1.5 and 2.0 Gyr. Afl sim- 
ulations evolve in a similar fashion: the initial gas distri- 
bution cools down slowly, loses pressure support and con- 
tracts in the vertical direction. After a few 100 Myr the 
central part of the disc is cold enough to become gravi- 
tational and thermally unstable and fragments into bound 
clouds. This process quickly proceeds to larger radii. Non- 
axissymmetric instabilities such as swing amplification aids 
the process everywhere, especially in the outer parts where 
the gas is only mildly unstable. The formation of bound 
clouds and elongated structures such as shearing filam ents 
is very similar to that found by e.g. iKim et aP (|2003l l and 
iKim Ostrikerl (|2007l ) for an unstable or marginally stable 
ISM. This clumpy structure is also visible in the evolution of 
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Figure 3. Logarithmic column density plots of the gas in the range Eg = lO"*^^ — 10^^ cm ^. Each panel shows a face-on 30 x 30 kpc"^ 
map centered on the disc. The associated edge-on map is 8 kpc in height. From top to bottom we see RUNl, RUN2 and RUNS at times, 
from left to right, t = 0.5, 1.0, 1.5 and 2.0 Gyr. 




Figure 4. Time evolution of the surface density (left) and rotational velocity (middle) for the gas component in RUNl. The contributions 
to the rotational velocity (solid line) at t = 1.0 Gyr (right) from the gaseous (dotted line), stellar (long-dashed line) and dark matter 
(dash line) components are in good agreement with observations of M33. 
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Figure 5. Time evolution of the mass of the gas and stellar com- 
ponent in a 30 kpc cube centered on the disc, as seen in Fig[3] 

the total gas surface density and rotational velocity in Fig.|4l 
The decrease of mean surface density is due to star forma- 
tion. Fig. U also shows the contribution to the rotational ve- 
locity at t = 1.0 Gyr from gas, stars and dark matter. We 
can clearly see that the initially gas and dark matter only 
system has evolved to a state in which the relative contribu- 
tion s and their m agnitudes agree well with M33 observations 
re.g lCorbellil[2003 ). The evolution of the total gas and stel- 
lar mass in RUNl and RUN3 are shown in Fig. [5] We note 
that a long-time evolution (t > 1.5 Gyr) of the galactic discs 
will force them to move away from a gas rich system such 
as M33, having ~ 30 % — 50 % of its baryonic mass in gas, 
approaching ^ 20% — 22% at t — 2 Gyr. Taking gas in- 
fall into account and using a more realistic star formation 
prescription could remedy this. 

We observe significant cloud-cloud and cloud-ISM in- 
teractions as the disc evolves. The clouds undergo both col- 
lisions leading to coalescence as well as tidal and long range 
interactions inducing torques into the gas. Shearing wavelets 
form out of the disc in between the cold clouds. These struc- 
tures interact with each other as well as the clouds for the 
entire simulati on period. The clum py ISM acts as an effec- 
tive viscosity (|Lin &: Pringld [ToStI ) forcing material to sink 
to the center whilst smaller clumps stay on more regular 
orbits at larger radii. This effect is stronger at early times 
when the typical cloud collision times-scale is short. In be- 
tween dense clouds and stretched filaments, the ISM also 
develops under-dense regions (E < lO^^cm"^) on scales of 
500 pc to several kpc. At later times, signatures of large scale 
spiral structure appear in the gaseous disc in which the cold 
clouds align. 

In what way do the simulations differ? RUNl and RUN2 
evolve in a very similar fashion. However, the higher resolu- 
tion in RUN2 means that further swing amplified instabili- 
ties can occur in the outer parts (see i — 0.5 Gyr in Fig. [3]). 
Apart from this, the overall morphology and statistics are 
in good agreement throughout the whole simulation time, 
indicating convergence. The feedback in RUN3 successfully 
ejects hot low density gas into the ISM as well as out of 
the disc plane. This process is very efficient at early times 
when the SFR/Area is high but calms down as the SFR self- 
regulates, see Sect. 13.3.11 for discussion. The feedback also 
alters the structure of the disc. As seen in Fig.El the late 



time spiral patterns are not as pronounced as in RUNl and 
RUN 2 and fewer low- mass clouds have survived. 



S.1.2 Star formation 

The left panel in Fig. [6] shows the SFR of RUNl, RUN2 
and RUN3. The main difference between the simulations is 
found in the most active star forming time, t ^ 0.2 — 0.5 Gyr, 
when the initial gravitational and thermal instabilities have 
formed dense clouds. The higher resolution in RUN2 allows 
for cloud formation in the less dense outer parts of the disc, 
leading to a higher SFR. The SNe feedback in RUN3 damp- 
ens star formation as explosions heat or disperse star form- 
ing clouds. However, after this star-bursting period the disc 
settles to a quiescent phase where all simulation approach a 
SFR - O.25M0yr-\ 

The stellar surface density is at all times well-fitted 
with an exponential function, X]*(r) ~ exp(— r/ro) with 
ro ~ 2.5 kpc, out to a truncation radius which grows with 
time, see right panel in Fig. (6] This may not come as a sur- 
prise as the initial gaseous disc is set up to be an exponential. 
At the end of the simulation time there are 1.3 x 10^ stars 
in RUNl, 8.8 x 10^ in RUN2 and 2.5 x 10^ in RUN3. 



3.2 Composition and state 

3. 2 A A multiphase ISM 

Theoretical models of the ISM fe.g. lMcKee k Ostrikel 19771 ) 
have a three phase structure consisting of a cold, warm and 
hot phase in pressure equilibrium where regulation is ob- 
tained through the balance of radiative cooling and super- 
novae heating. More updated models separate the phases 
further based on their ionization state. In this work we re- 
fer to the phases as cold (T < 10^ K), warm (10^ K < T < 
10^ K) and hot (10^ K > T). A realistic ISM is very compli- 
cated which can be seen in the volume-weighted phase dia- 
grams of RUNl and RUN3 in Fig.fT] Both runs show a wide 
range of temperatures and densities. RUN3 clearly displays 
distinct cold, warm and hot phases aligning to an isobaric 
strip, i.e. P ^ pT ^ constant. Around this region we observe 
a large spread in both temperature and density. This analy- 
sis is approximately valid for RUNl even though the warm 
and hot gas are smeared over less distinct phase regions and 
sits at slightly lower densities compared to RUN3. The cold 
phase is almost identical in the two models. The existence of 
a hot tenuous phase in RUN3 is due to supernovae heating, 
but why do we find hot gas in RUNl? As the initial circular 
velocity is set for the whole computational domain, the low 
density ambient gas at T = 10^ K experiences a mild shock 
heating and settles into pressure equilibrium with the denser 
galactic disc. Some of the hot gas can also be found in the 
cloud-induced shocks in the disc. 

Fig. [8] shows the volume fraction occupied by gas at dif- 
ferent temperatures. We also indicate the total percentage 
of gas in different temperature regions. The expected two- 
phase structure in RUNl is evident. We find a clearly peaked 
cold phase with a transition into a warm phase in between 
1000 K and 10 000 K, peaking at 6 000 - 7 000K which is 
the thermally unstable regime. The origin of the warm phase 
is shock-heating. RUN3 shows the same cold phase but the 
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Figure 6. (Left) The star formation rate over time for RUNl, RUN2 and RUN3. The higher resolution in RUN2 allows for more star 
formation in the outer disc, while the feedback of RUN3 lowers the efficiency. At late times the SFRs show little difference. (Right) 
Evolution of the stellar surface density. The density profiles are at all times well fitted with an exponential function. The red line is for 
ro = 2.5 kpc. 
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Figure 7. Phase diagrams for RUNl (left) and RUN3 (right) at t = 2.0 Gyr. The solid straight line indicate the isobar. 




Figure 8. Volume fraction for RUNl (left) and RUNS (right) at t = 2.0 Gyr. While both models clearly show a cold and warm gas 
phase, the hot phase is only present in RUN3. 
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Figure 9. Most unstable wave lengths for the whole disc around 
time of fragmentation i.e. t ~ 100 Myr. A large part of the disc 
is unstable at ~ 100 pc — 2 kpc, scales that will collapse to the 
initial cloud distribution. 



warm phase now strongly peaks at 10 000 K, just at the max- 
imum peak of the cooUng function. A hot gas phase is clearly 
present even though very little gas exists above 10^ K. As in 
RUNl, the warm phase dominates the gas volume. 

It is desirable to approximately reproduce a mass distri- 
bution of molecular, cold atomic, warm atomic, warm ion- 
ized and hot ionized gas that agrees with observations (see 
e.g. iFerrierd (|200lh for the Galactic inventory). However, 
even among the local group spirals there can be significant 
differences between the phase-distributions. For example, 
M31 has - 40% of its gas in cold HI, Milky Way - 25% 
and M33 only - 15% (Dickev k Br inks| 1 19931). These dif- 
ferences could be due to the variation in baryon to dark 
matter fraction as we move down the Hubble sequence. As 
we show later, the formation of cold clouds is particularly 
sensitive to the gaseous disc mass. However, it is still in- 
structive to compare our ph ase value s of RU Nl and RUN3 
at t = 1.5 Gyr to those of IFerrierd (|200lh for the Milky 
Way. Roughly 50 % of the Milky Way gas is in molecular 
and cold atomic (50 — 100 K) clouds. Our simulations only 
allow for cooling down to 300 K and can hence not dis- 
criminate between the coldest gas phases. By labeling all 
dense gas of T < 350 K as a joint cold cloud phase we 
find that - 55(47) % of the total gas mass in RUNl (3) is 
cold. The warm neutral gas phase (10^ K < T < 10"^ K) has 
^ 11(16) % while the total neutral mass fraction of the gas 
outside of clouds (350 K < T < 10^ K) is 39(45) % respec- 
tively. The former value is lower than the Milky Way value 
40 %) which could be due to the fact that our initial con- 
ditions are more suitable for comparison with Sc galaxies. 
Also, including a homogenous UV background field in the 
simulations would heat the dif fuse HI gas wh ich seems to be 
the case in similar studies fe.g. lBottemall2003l l . Furthermore, 
the observed mass of warm phase in the Galaxy is derived 
from the observed HI velocity dispersion of 6 — 9 kms^ 
unde r the assumption of only thermal broadening (|Ferrier3 
l200lh . A turbulent component can allow for the existence of 
colder gas yet retaining the velocity dispersion values. We 
will explore this notion further in Sect. l3.3Tl 
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Figure 10. Cumulative mass spectrum of individual "molecular" 
clouds (n > lOOcm-3) in RUNS. 



3.2.2 Disc stability 

To understand the relevance of gravitational instabil- 
ities in the simulat ed mult ip hase discs, we use the 
Toomre parameter (Flbomrd [l964[ ) defined, for gas 
(|Goldreich Lvnden-Belllll965al ). as 



Qg 



(8) 



where Cs is the sound speed of the gas. Since the gas has 
turbulent motions it is appropriate to use the effective dis- 
persion (Jeff = Cs + cTiD, where ctid is the average of the 
full three-dimensional velocity dispersion. The Toomre pa- 
rameter is valid for local axisymmetric perturbations of two- 
dimensional discs, where Qg < Qc = 1 implies instability. 
However, Qg has been shown to characterize the response 
of discs to general gravitational instabilities. A finite disc 
thickness weakens the surface gravity and lowers the criti- 
cal value where the disc u ndergoes instability. For example, 
iGoldreich k Lvnden-Bell ([^65a) showed that Qc = 0.676 
for a single-component thick disc. In addition, the onset for 
non- axis symmetry occurs at higher values of Qg, both for 2D 
(Qg ~ 1.7) and 3D discs. Extended stability analysis tak- 
ing thickness and multiple compo nents (collisional and /o r 
collisionless has developed b y e.g. IJog Solom"^ (|l984l l: 
iRomeol (|l992l ): iRafiko^ (|200lh . In this section we mainly fo- 
cus on the more unstable gas component as it will be the 
main driver of turbulence compared to the stellar compo- 
nent which shows a higher degree of stability at all times 
{Q. > 2). 

We start by investigating the early time evolution when 
the initial cloud population forms. To do this we need to 
quantify the most unstable length scales. The dispersion re- 
lation for axisymmetric disturbances 



2 2 
UJ — 



(9) 



where uo is the growth rate and k is the wavenumber of the 
perturbation. Instability demands that <0 and the most 
unstable mode is simply the minima of Eq.[9l i.e. 



An 



(10) 



We calculate all components the dispersion relation, and 
Qg, on a polar grid with subregions 30kpc/Ai? = 40 and 
27r/A^ = 60. Fig. [9] shows the minimum of the dispersion 
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Figure 11. (Left) Volume distribution of Qg values for the whole disc in RUNl. The low Qg values at early times indicates the formation 
and existence of clouds. Later times shows an equilibrium distribution that changes little with time. (Right) Mass distribution of Qg 
values. Most of the disc is distributed over 1 < Qg < 4. 



relation and their related growth factors at t = 100 Myr. 
The bimodal distribution of points might seem odd but is 
merely a reflection of the initial Q(r) (see Fig. [2]), where the 
same value of Qg exists at different radii (and different den- 
sities) and hence show the same at different Amin- We 
note that scales of 100 pc < Amin < 2 kpc are in the un- 
stable regime, where the smaller scales show larger negative 
values of co^. It is reassuring that smaller scales remain sta- 
ble due to the imposed Jeans capturing EOS discussed in 
Sect. 12.41 These gravitational instabilities set the initial con- 
ditions for the clouds. The evolution of the cumulative mass 
spectrum, for "molecular" gas (n > 100 cm~^), in RUNl 
is shown in Fig. 1101 We have calculated the mass spectrum 
by simply discerning individual pieces of high-density gas 
in the disc. This method is crude and occasionally overes- 
timates the mass of clouds in the central parts of the discs 
where the gas density is high and crowding artificially iden- 
tifies several clouds as one. Disregarding this, we note that 
the spectrum around t ^ 1.5 Gyr occupies sim ilar values 
as that of local group spirals (|Blitz et al.ll2007h , where the 
most massive clouds are ^ 10^ Mq. The small mass trun- 
cation is due to limited resolution. As the simulations lack 
important small scale physics, e.g. MHD, radiative trans- 
fer, cosmic rays etc., the cloud population is long-lived and 
only reflects a true ISM in a statistical sense. This notion 
should not be a problem for the source of turbulent veloc- 
ity dispersions that, as shown in Sect. 13. 3[ is due to large 
scale gravitational drags that most probably is independent 
of the small scale gas state close to or inside of the cloud 
complexes. 

We now turn to the subsequent evolution. Fig.fTTIshows 
the time evolution of the distribution of Qg values for the 
whole disc. The left panel shows the volume fraction that 
different values of Qg occupy while the right panel treats 
the mass fraction. This flgure illustrates the complexity of 
the simulated discs and why azimuthally averaged Qg(r) 
can be misleading. Initially, the disc shows a low spread of 
Qg around a value of a few. As the disc cools down and 
undergoes gravitational instability (after t ^ 0.1 Gyr) this 
simple picture changes. At t = 0.25 Gyr, the disc has un- 
dergone fragmentation and the distribution is confined to 
0.2 < Qg < 1. The peak of the distribution, and the dis- 
persion, gets larger with time. Part of this owes to star 



formation that acts to lower Eg. For t > 1.0 Gyr the disc 
evolves into what appears to be an equilibrium state, span- 
ning a large range in Qg-values (0.5 ^ Qg ^ 10^). This 
co-existence of Qg-value in a patchy galactic disc is in agree- 
ment with the analysis of W ada et al.. (.2002i ) for their two- 
dimensional models. The mass fraction distribution follows 
a similar evolution, approximately reaching an equilibrium 
state after t > 1.0 Gyr. However, a significant part of the 
Qg distribution at late times now populates the unstable 
or marginally stable values. At 1.5 Gyr, most of the disc is 
distributed around 1 < Qg < 4. Regardless of the exact dis- 
tribution, the dominating existence (by mass) of unstable 
and marginally stable regions of the disc is of great impor- 
tance for generating a global gravitoturbulent state which 
we will return to in Sect. 13.3.21 Without the onset of grav- 
itational instabilities, the gas would approximately stay on 
circular orbits. 

At late times the mass in the disc is dominated by the 
stellar component. To get an understanding of its infiuence 
on stability one can use an approximat e stability parame- 
ter (jBertin Romeolll988l : lRomeolll994l ) which in standard 
regimes is of the form 



Q~Q* 1 



(11) 



At late times, the stellar component is in the range 1.5 < 
Q* < 5 in the star forming region which together with the 
values shown in Fig. [11] assures us that the multicomponent 
disc will never be completely stable, at least locally. 

3.3 The turbulent ISM 

Having characterized the multiphase ISM in terms of phases 
and stability, we can now properly address the main topic 
of this work, the HI velocity dispersions. 



3.3.1 Velocity dispersions 

The observational tradition is to model HI profiles using 
one or multiple Gaussians where the fiux is a function of the 
velocity v as 



f(v) 



aV2 



exp 



■ vo) 



(12) 
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Figure 13. Observed velocity dispersion (blue dashed line) of 
the face on galaxy NGC 1058 compared to the effective vertical 
dispersion (black solid line) of RUNl at t = 1.5 Gyr. Grey dots 
indicate locally measured dispersions, see text. 



where vo is associated with the peak flux and a is the actual 
velocity dispersion. Broadening of spectral lines is mainly 
due to thermal and Doppler/turbulent effects. We will dis- 
cuss the thermal effects in terms of the thermal velocity 
vt — ^JRT / fi (i.e. the isothermal sound speed of the gas), 
where R is the gas constant and /i is the molecular weight. 
Random bulk motion of the gas is quantified in terms of 
its turbulent velocity dispersion at. We calculate the net 
observable dispersion by adding the turbulent and thermal 
contribution in quadrature, i.e. deff — + (Jt • 

In the following analysis, we characterize only the gas 



that would be observed as HI and not the dense clouds 
that will consist of mainly molecular gas. We therefore use 
the criteria p < 10 cm~^ (star formation threshold) and 
800K < T < 10 000 K. This choice is suitable as it is more 
Ukely to exist outside of the denser spiral arms as well as 
in the outer regions of the disc, which is where observations 
lack an explanation. The velocity dispersion is calculated 
by randomly sampling the galactic discs using synthetic ob- 
servational patches (5 000 patches were used for the data 
described here) of size ^ 700 pc. We choose this size as this 
is the stated scale below which bulk motions are expecte d to 
be re sponsible for the observed dispersions (Petric & R upenI 
l2007l ). In each patch we calculate both the mass weighted 
turbulent velocity dispersions at and the mass weighted 
mean thermal velocity. Weighting by mass is well motivated 
as HI emission is strongly correlated with the local density. 

The panels in Fig. 1121 shows a time evolution of the ra- 
dial behaviour of the velocity dispersion for RUNl, where 
Gz is the vertical dispersion component, Gr the radial and 
(T(f) the angular. (Tz show typical values of 15kms~^ in 
the center and declines to ^ 3 — 5 kms~^ at large radii. The 
velocity dispersion is clearly anisotropic as dr > cr^ > Gz at 
all times. It is interesting that the ratio of the dispersions 
roughly follow the epicyclic predictions for a collision-less 
syste m, i . e. Gr — 2Q.G^/Hi. A similar result was found by 
BottenS (l2003l). The planar dispersion Gxy^ i.e. the RMS 
value of the radial and angular dispersions, is a factor of 
2 larger than the vertical dispersion at all times and 
radii. The thermal component of the gas lies in the range 
3 — 5kms~^ in agreement with a warm gas component 
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(T ~ 1000 — 2000 K). The planar velocity dispersion is super- 
sonic or transsonic at all times and radii while the vertical 
dispersion is generally transsonic, turning sub-sonic at large 
radii. This means that the thermal component becomes as 
important as the turbulent at large radii for the total observ- 
able velocity dispersion. By considering a minimal observ- 
able (fJeff) for the z-component, we clearly find an agreement 
with the observed HI dispersions values described in Sect.[T] 
(i.e. (Jeff ^ 12 — 15 kms~^ in the inner parts declining to 
~ 4 — 6 kms~^ in the outer). Any inclination would boost 
these values due to the cr-anisotropy. The same analysis has 
been performed on the higher resolution simulation RUN2 
with no significant difference in the results. 

We now directly compare our si mulations to the HI 
data of the spiral ga laxy NGC 1058 (IPickev et al.1 I199QI : 
IPetric RuDenll2007l ). NGC 1058 is a suitable object for 
comparisons, as it is comparable in size, surface density 
and peak rotational velocity (derived to be ^ 150kms~^) 
to our simulated disc. The galaxy also has a low star for- 
matio n rate, SFR~ 3.5 x lO~^M0yr~^ dPerguson et al.1 
11998V which places it in the flat part of Fig.[T] Furthermore, 
by being an almost perfectly face-on galaxy (inclination of 
4— 11 °), we can disentangle the vertical component from the 
planar. In Fig. [13] we compare cr2,eff of RUNl at t = 1.5 Gyr 
with the observational data of NGC 1058. Our simulation 
not only reproduces the magnitude of the velocity dispersion 
but also the declining radial shape. 

The spatial distribution of the vertical velocity disper- 
sions in NGC 1058 is very patchy with several p eaks of 
a > lOkms" (see Fig. 5 of lPetric Ruuenl (|2007l )). This 
observation is reproduced by our simulations, as shown in 
Fig. [14] where we plot the contours of creff,^ in RUNl and 
RUN3 at t = 2.0 Gyr as well as the corresponding density 
fields. The high density gas is distributed in a flocculant 
spiral structure, reminiscent of the HI observation of M33 
(jPeul van der HulstI Il987l : lEngardola et al.ll2003l ). Ana- 
lyzing the simulation at a late epoch is preferred as the mass 
spectrum of dense clouds has evolved to a rather realistic 
state, see Fig. [10] The strongest peaks {(Jef£,z > lOkms""*^) 
are associated with dense clouds while mildly turbulent re- 
gions (cyan levels at creff,^ 6 — 8kms~^) often exist in 
inter- cloud/ arm regions of strong shear. Regions of large ve- 
locity dispersion related to clouds also extend several kpc 
away from their radius of influence. We note that RUN3 even 
at late times displays a few kms~^ larger velocity disper- 
sions in diffuse regions, probably due the more prominent 
warm gas phase. 

3.3.2 What is the driver of ISM turbulence? 

The drivers of the turbulent component of the velocity dis- 
persion are gravity and shear. In Sect. 13.2^2] we found that 
the galactic discs have, by mass, a wide spectrum of Qg val- 
ues where a significant part sits at local marginal stability 
for a finite thicknes s disc. The 2D shearing box simulations 
by iKim &: Ostrikerl (|2007h showed that a marginally stable 
gas discs at Qg ^ 1.2 can generate velocity dispersions of 
the order of the local sound speed, decreasing for larger Qg- 
values (see their Fig. 12). The origin of turbulence was here 
attributed to swing amplification. Note that the 3D struc- 
ture of our discs necessitates values ~ 25 % lower for an 
equivalent stability. A full turbulent outcome of more un- 



stable discs (Qg < 1.2) was not studied as the velocity field 
would then only be a response to very strong density inho- 
mogeneities. Local shearing boxes are useful for understand- 
ing the mechanism that drives turbulence at specific values 
of Qg. As our simulations show a wide range of Qg values 
we have the combined spectrum of swing amplified turbu- 
lence across the whole disc for gas that locally beh a ves in 
accordance with the simulations of iKim Ostrikerl (|2007f ) 
for Qg > 1.2. In a statistical sense there will always be re- 
gions with a Qg-value low enough to tap large velocity dis- 
persions from the swing mechanism which is confirmed in 
Fig. [14] where intermediate values of aef£,z is associated with 
waves. To qua ntify this it is useful to use the X-parameter 
(|Toomrelll98"ll ) defined as 

X = (13) 
m 

where /ccrit — A€^/27rG'5], S = Eg + arid_m is the num- 
ber of arms. It has been shown by Ijod (|l992l ) that swing- 
amplification is very effective in the gas component in mul- 
ticomponent discs. Even when both the gas and stars sep- 
arately are stable (Qg = Q* = 2), their gravitational cou- 
pling can amplify waves in he gas for values of X not much 
larger than unity. In one component, marginal stability and 
1 < X < 3 can be considered sufficient to assure amplifica- 
tion (|Toomrelll9"8ll ). Fig.l 151 shows the radial dependance of 
the X-parameter for m = 2, 4 and 8 at t = 2.0 Gyr. We see 
that amplification is efficient for m > 4 which confirms the 
high-order flocculant spiral structure in Fig. 1141 

For small values of Qg, where the gas locally has under- 
gone full non-linear gravitational instability, the situation is 
different. The cold phase dominates the gas mass, even at 
early times and is therefore locally the most important grav- 
itational source. Direct cloud merging and tidal interactions 
stirs the inter-cloud medium both radially and vertically. 
Apart from stirring the gas, the clouds also dissipate energy 
thermally in shocks which regulates the warm phase of the 
ISM, forming the 4 — 5kms~^ thermal components of 
(Jeff. We observe cloud formation, merging, scattering and 
reformation during the whole simulation time. Formation in 
a shearing environment causes dense structures that are not 
tightly bound to the actual clouds to stretch into waves and 
filaments. This triggering of wave-like perturbations, and its 
associated irregular velocity field in the ISM, is a key ro le of 
the c louds which was realized already by Julian Toomrd 
(| 19661 ). As for the marginally stable gas surrounding the gas, 
the leading waves swing and amplify, inducing gravitational 
torques in the gas and h ence incr easing the local velo city dis- 
persion fe.g. lKim et HI [20021 : IKim Ostrikerll2007l ). Fig.[T6] 
shows a typical patch of the disc in RUNl at t = 1 Gyr, con- 
firming this notion. We note that this processes is analogous 
to the energy extraction from background shear at a rate 
TR(f)dQ / d\n where the Tr^ stress tensor includes the con- 
tribution from Reynolds and Newtonia n stresses, to induce 
local velocity dispersion as outlined by ISellwood k BalbusI 
([l99i). 

But how can we quantify the impact of the cloud mo- 
tions? Let us assume that the motions of the cloud ensem- 
ble are representative of the turbulent ISM. The swing am- 
plifier might play a very fundamental role for turbulence 
but as the clouds effectively trace the large scale waves and 
constitute the majority of the mass, the assumption that 
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Figure 14. Distribution of vertical velocity dispersions in RUNl {top left) and RUNS (top right) at t = 2.0 Gyr calculated for HI only 
(see text). The plotted region is 30kpc across. The contour levels are in kms"-*^ in steps of Ikms"-*^. Black is used between 3 and 5, 
cyan for 6 to 8, red for 9 to 11 and green for 12 to 14 kms"-*^. The bottom panels show contours of the gas with n > 0.2 cm~^ of the 
same regions but slightly zoomed out (40kpc across). 
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Figure 16. Density plot of a 24 x lOkpc^ region centered over {x,y,z} = {0, — 5.0,0} kpc in RUNl at t = 1.0 Gyr. The colour map is 
here chosen to enhance the visual appearance of the clouds and filaments. The galactic rotation is here clock-wise. Filamentary structures 
is always associated with the clouds. Also, all clouds excite waves, many of them leading which will swing into trailing ones. 




^ 4 6 13 IQ 1^ 

r fkpc] 

Figure 17. Planar velocity dispersion of RUNl at t = 2.0 Gyr 
(black solid line) compared to the relation derived by Gammie et 
al. (1991) (blue dashed line) for cloud scattering. The epicyclic 
frequency is obtained from the simulation and the cloud mass is 
chosen to be 3.5 x 10^ M©. 



turbulence is associated with cloud motions is a fairly good 
approximation. Cloud-cloud interaction can be modelled as 
gra vitational scatter i ng; an d has been studied analytica lly by 
e.g. Ijog Ostrikerl (|l988l l and lGammie et al.1 (Il99lh. The 
semi-a nalytical perturbation theory model of iGammie et al.] 
(|l99lh predicts a planar velocity dispersion 



,1/3 



(14) 



where Md is a typical mass of a cloud. This relation is de- 
rived for a two-dimensional, two body encounter on radi- 
ally separated orbits in a shearing disc. However, as these 
clouds are the main local perturbers by mass we can as- 
sume that any diffuse HI gas will approximately be dic- 
tated by the cloud ensemble velocities. In Fig.[T7l we plot 
axy for RUNl at t = 2.0 Gyr against Eg. [Ml using a cloud 
mass Mci ~ 3.5 x 10^ M©, and /^(r) of the gas in the sim- 
ulation. We stress that Md is in the high -end of a typi- 
cal GMC mass spectrum (Blitz et al. 2007). As the clouds 
in our simulations are submerged in massive HI envelopes, 
the largest clouds are closer in mass to that of GMC com- 
plexes, CM As (Giant Molecul ar Associations) or super- 
clouds (iRosolowskv et all l2007l ) . A more realistic analysis 
should include the full spectrum of cloud masses and their 



radial distribution but even this simple analysis renders a 
good agreement with the measured dispersions. The weak 
dependence on k, can explain why most non star-bursting 
galaxies seem to plateau at a velocity dispersion between 7 
and llkms"^ (see Fig.[T]and discussion in Sect.[T]). The ro- 
tational velocity varies from ^ lOOkms"^ to ^ 300kms~^ 
for most spirals. By assuming a flat rotation curve, k oc Vc, 
the actual change in cloud velocity dispersion between the 
two Umits is, assuming an invariant cloud spectrum, only by 
a factor of 3^/^ 44%). 

To conclude, the full picture of gravity driven turbu- 
lence in the ISM is based on the existence of dense clouds and 
filamentary structures all adding to the turbulence budget. 
Marginally stable gas and cloud-induced filaments generate 
turbulence through gravitational torques from the swing- 
amplifier. Clouds also scatter gravitationally, perturbing the 
local velocity field even further as well as shock-heating the 
ISM. The source of turbulence in both cases is self- gravity 
coupled with shear that in turn converts ordered circular 
motion of the gas to random velocities, hence tapping rota- 
tional energy from the disc. In a s ense, self- gravity can here 
be regarded as a form of viscosity (|jog Ostrikerlfliiih . We 
argue that the mechanisms described in this section serve as 
a baseline level of turbulence for galaxies where any excess 
observed velocity dispersion is caused by additional sources 
such as supernovae activity or magnetic fields coupled with 
shear. Finally, we note that many of the disc characteristics, 
e.g. a and Qg, show signs of reaching a statistical equilibrium 
state at late times in the simulations. This is an indication 
that there is a constant supply rate of energy to the system, 
here coming from galactic rotation, maintaining a constant 
level of turbulence. However, on close inspection these quan- 
tities are slowly decaying (e.g. the mean velocity dispersion 
is slowly decreasing, see Fig. ll2|) due to star formation de- 
pleting the disc of gas. Inclusion of realistic gas accretion 
would affect the temporal evolution of the latter quantities. 



3. 3. 3 Varying the baryon fraction, shearing motions and 
cooling floor. 

In order to explore the effect of the gas density in the initial 
disc, we keep all parameters fixed apart from the disc mass 
and carry out an additional simulation (RUNG) in which 
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Figure 18. Cumulative mass spectrum of clouds in RUNO 
(dashed line), RUNl (black) and RUN6 (red) at t = 0.75 Gyr. 
The smaller region of instability in the less dense disc in RUNG 
brings down the total number of clouds and the mass offset is 
related to a decrease of the largest wavelength to be stabilized by 
shear, see text for discussion. 



Figure 20. Ratios of epicyclic frequencies for RUNl, RUNG and 
RUNS. The RUNl/RUNG ratio is above unity in the central parts 
owing to excess mass transport to the center in RUNl and hence 
an increase in central shear. The overall increase in RUNS, where 
the central parts of the disc are more affected, is expected from 
the change in dark matter halo concentration. 




r [kpc] 

Figure 19. Ratios of planar velocity dispersions for RUNl, RUNG 
and RUNS. The RUNl/RUNG ratio illustrates the impact of the 
surface density of the gas, where the more massive disc in RUNl 
has ~ 3 times larger velocity dispersion. The increased shear in 
RUNS increases the velocity dispersion when compared to RUNl, 
but only significantly in the central parts of the disc. 
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Figure 21. Ratios of planar velocity dispersions for RUNl, 
RUN4, RUN5 and RUNT showing the effect of a temperature 
floor for the gas cooling. The velocity dispersion in RUN5 and 
RUNT are lower by a factor of 3-4 compared to RUNl due to 
the inability to undergo gravitational instability. The RUN4 disc 
is allowed to cool enough to capture the largest clouds forming, 
hence having similar velocity dispersions as RUNl. 



the gas mass is one third of our fiducial value. This changes 
the surface density profile and thus the strength of self- 
gravity. This disc is extremely light and is unstable only 
within ^ 8 kpc, leading to an abundance and mass spec- 
trum of cold gas clouds that is significantly smaller than in 
RUNO or RUNl as seen in Fig. [18] Clouds are here defined 
as isolated clumps of gas satisfying n > 100 cm~^. The ver- 
tical shift, when comparing RUNO and RUNG, is mostly due 
to the larger area of instability in the latter simulation and 
the RUNO and RUNl off- set is simply due to star- format ion 
acting to deplete the clouds of high-density gas. We can 
quantify the observed mass-shift using the arguments pre- 
sented by e.g. .Escala Larson (2008) where the high-mass 
end of the spectrum is set by the largest modes not to be 
stabilized by shear, i.e. the maximum cloud mass will be set 
by 



4^4 



(15) 



A factor of 3 decrease in Sgas therefor leads to M™^"" be- 
ing 27 times smaller which is in excellent agreement with 
the mass spectrum shift in Fig. [18] The most massive clouds 
in RUNG are ^ lO^M©, leading to weaker cloud-cloud en- 
counters and swing amplified waves, directly lowering the 
accelerations imparted on nearby gas parcels. This leads to 
a lower overall velocity dispersions by a factor proportional 
to M^^ - S (from Eq. [Hand Eq.^. An off-set by a fac- 
tor of ^ 3 is confirmed in Fig. [19] which shows the ratios of 
o'xyir) of the different simulations. 

To isolate the effect of the shearing motions we keep 
all parameters the same as our fiducial RUNl, apart from 
the rotation curve set by the dark matter halo which is con- 
structed to be much flatter (RUNS). This is achieved by 
increasing the concentration parameter to c = 40 whilst 
maintaining the peak circular velocity to be similar to the 
rest of our simulations. The swing instability and turbu- 
lence induced by shearing motions should be stronger since 
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the epicyclic frequency, defined as 



(16) 



wfiere Q = Vc/r^ increases. For a Vc 
K, ~ and for Vc '^constant, k 



yjr this mean that 
1/r. The former Vc 
roughly describes RUNl and the latter RUNS. Fig.[20lshows 
K(r) for RUNl, RUNG and RUNS at t = 0.75 Gyr. RUNl 
and RUNG show similar values apart from the inner part 
of the disc where turbulent viscosity in RUNl has dragged 
in more mass compared to RUNG, rendering a more active 
region. RUNS shows a steeper behavior, having a At up to 2.5 
times larger in the very center of the disc down to a ratio of 
unity at r ^ Skpc. The larger shear causes a larger (Jxyir) 
compared to RUNl, as seen in Fig. 1191 The effect is strong 
in the central parts of the disc (1.5 — 2 times larger) and for 
r > 2 kpc the ratio is closer to unity. 

It seems evident that gravitational instability is the im- 
portant driver of velocity dispersions. To explore this fur- 
ther, we artificially truncate the ability to cool gas at differ- 
ent temperatures, hence effectively setting a floor for Qg . 
In RUN4, RUN5 and RUNT we introduce a cooling floor for 
the gas at 10^, 10"^ and 5 000 K respectively. In Fig. [211 we 
compare the ratios of Gxy(r) at t = 0.5 Gyr and it is evident 
that RUN5 and RUNT only show small turbulent motions 
while RUN4 essentially is as turbulent as the fiducial RUNl. 
The physical parameters are the same in all simulations ex- 
cept for the temperature of the gas. As Qg a '-^ \/T it 
is straightforward to see that RUN5 and RUNT, based on 
the initial Qg ~ 2 — 3 in Fig.[2l never acquires Qg < O.GTG. 
Both simulations only show a weak development of spiral 
structure. We note that it is plausible that a stellar compo- 
nent wo uld increase the turbule nt motions of the discs, as 
found bv lKim h Ostriked (|200TI ). Also, the effective observ- 
able cr^^eff in RUN5 and RUNT is dominated by the ther- 
mal component as is at large radii close the values found 
in RUNl. The cooling in RUN4 can lower Qg by a factor 
~ 3.2, hence bringing Qg < O.GTG while being assisted by 
non-axissymmetric instabilities. While the smallest scales 
to be unstable differ in RUNl and RUN4 due to different 
pressure support, the larger unstable wavelengths are the 
same, only limited by the same amount of shear. The pres- 
ence of the larger, and more dynamically important, clouds 
and gravitational instabilities makes the velocity dispersion 
ratio closer to unity. 

We conclude that the mass, and hence the surface den- 
sity, of the disc has a significant impact on the generated 
turbulent velocity field and we roughly find that, provided 
the disc is gravitationally unstable, a ^Y,. We associate this 
to the weaker gravitational instabilities present in the disc, 
which can be seen from the cloud spectrum. The shear of 
the disc also affects the magnitude of the velocity disper- 
sions but in a weaker fashion. This less strong effect of shear 
might originate from the ^ a^^^"^ -relation fEa. ll4|) for cloud 
velocity dispersion. 



S.S.Ji. Effect of supernovae feedback 

Self-gravity driven turbulence may be important for galax- 
ies with a low SFR/Area but cannot be the dominant driver 
behind the large velocity dispersions correlated with high 
star formation rates. Observations suggests that galaxies 
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Figure 22. Effect of star formation rate on the observed vertical 
velocity dispersion in RUNS. The different lines indicate the val- 
ues at different times and therefore also for different SFRs. There 
is a clear trend that a lower stellar activity lowers the measured 
dispersion, approaching the baseline observed dispersion given by 
RUNl at t = 2.0 Gyr (thick red dashed hue). 
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Figure 23. Effect of star formation rate on the observed vertical 
velocity dispersion in RUN3. The different lines indicate the val- 
ues at different times and therefore also for different SFRs. There 
is a clear trend that a lower stellar activity lowers the measured 
dispersion, approaching the baseline observed dispersion given by 
RUNl at t = 2.0 Gyr (thick red dashed hue). 



with a SFR/Area > few x lQ~^M(:)yr~^kpc~^ show ve- 
locity dispersions of several lOkms"^, see Fig.fTI [Pib et al.l 
(|2006l l showed that strong SN feedback could explain the 
transition into this range but were unable to explain the 
other end of the spectrum. It is plausible that the reason 
for this stems from their local shearing box approximation 
that does not take the full disc dynamics into account. Fur- 
thermore, IPib et al.l (|2006l l demonstrated that supernovae 
driven turbulence is sensitive to poorly known parameters 
such as efficiency, mass loading, timing etc. Due to com- 
putational cost, we can only study one set of parameters. 
However, as the star formation rate decreases with time (see 
Fig. [6]) we are able to study its correlation with velocity dis- 
persion. In Fig. [22] we plot the effective vertical dispersions 
for RUNS at different times and hence different SFRs. The 
general amplitude of the dispersion declines with SFR and 
after t = 1.5 Gyr (SFR ^ 0.74 Moyr"^) there is little dis- 
crepancy between RUNl and RUNS suggesting that the ef- 
fect of SN feedback has saturated. As seen in Fig[8l more 
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warm gas (close to ~ 10 000 K) exists in RUNS explaining 
the ^ lkms~^ off-set between RUNl and RUNS at large 
radii. This difference is also seen by inspection of the cyan 
contours in Fig. 1141 

The data shown in Fig.[T] can be reproduced by av- 
eraging the velocity dispersion and SFR over a suitable 
area, he nce obtaini n g a c haracteristic velocity dispersion 
Cchar (r). iDib et al.1 (I2OO6') used the area A — 7r(Sro) , 
where ro is the scale radius of the stellar disc. Fig. [23] 
shows the outcome of this procedure and we clearly de- 
tect a supernovae saturation to occur at a SFR/ Area of 
1-2x10"^ Mo yr"^ kpc"^ where (Jchar for RUNl and RUNS 
coincides at t = 2.0 Gyr. This transition from supernovae 
to self-gravity induced turbulence can explain why the ve- 
locity dispersion of NGC 1058 presented in Sect. lS.S.T] is in 
good agreement with our simulated disc. NGC 1058 has a 
derived SFR ^ S.5 x 10"^ Moyr"^ (IPetric k R upen 2007) 
which sets the SFR/Area weh below lO"^M0yr~^ kpc"^ and 
hence into the regime where self-gravity induced turbulence 
can explain the observations. 



4 CONCLUSIONS AND DISCUSSION 

Three-dimensional, high-resolution hydrodynamical simula- 
tions using realistic modelling of star formation and evolu- 
tion, show that a turbulent ISM naturally develops due to 
the coupling between gravitational instability and shearing 
motions. A multiphase medium develops in which cold dense 
clouds and filaments co-exists with a diffuse warm gas. When 
supernovae feedback is implemented, a hot phase is present. 
The marginally stable gas undergoes swing-amplification 
which both acts to amplify the local density as well as in- 
ducing gravitational torques. Cold and dense clouds undergo 
gravitational scattering, merging and tidal encounters. They 
also induce waves and filaments in the more diffuse gas which 
pumps energy into the turbulent process. The former mech- 
anism stirs the gas even further and we note that the veloc- 
ity dispersion of the clouds is a fairly good tracer of the HI 
velocity dispersion. 

We summarize our main conclusions here: 

• Gravitational instabilities in galactic discs leads to a 
population of massive cold clouds that undergo mutual grav- 
itational interactions and merging. This cloud-cloud harass- 
ment process strips material and stirs the ISM. Both cloud 
interaction and the global non-axisymmetric instability of 
the disc create non-circular motions from initial ordered ro- 
tation. Waves and filaments are generated in the ISM which 
in turn swing- amplifies to generate further turbulent mo- 
tions. 

• Below a star- format ion rate per unit area of 
lO~^M0yr~^ we find that gravity alone can provide the en- 
ergy source for maintaining the observed level of turbulence 
in the ISM of galaxies. The turbulent velocities in our MSS 
model galaxy have a mean value of ^ lOkms"^. By calcu- 
lating an observable HI velocity dispersion, i.e. the contribu- 
tion from both the turbulent and thermal components, we 
show that both the magnitude and radial profile is in good 
agreement of high- r esolution HI surveys o f e.g. NGC 1058 
(iDickev et al.lll990l : [Petric RuDenll2007l ). In addition, we 
reproduce the observed patchy velocity dispersion map. 



• Once the star- formation rate exceeds this value, su- 
pernovae feedback becomes the dominant driver of turbu- 
lence and the velocity dispersion increases with the star- 
formation^^ e^^^is agrees well with the general trend found 
bv lDib et al.1 (|2006i) . 

• Lowering the initial gas density weakens the strength 
of gravitational instability and lowers the resulting cloud 
mass spectrum, which in turn leads to a lower disc velocity 
dispersion by a factor oc M^Y^, as expected from a model in 
which self-gravity generates significant turbulent motions. 

• A direct prediction of this scenario is that galaxies with 
lower gas fractions at a fixed halo mass should have lower ve- 
locity dispersions and different mass fractions in cold, warm 
and hot phases. Although, detecting the dependence on sur- 
face density is complicated by the fact that lower mass galax - 
ies have a higher gas fraction in their discs ( McGaughl l2005l ) . 
In addition, the reaction in low-mass systems to mild stellar 
activity has not been tested in this work and the outcoming 
HI velocity dispersion might conspire to render the plateau 
in Fig.m 

It is important to note that these results do not rule 
out the importance of other contributing mechanisms such 
as supernova feedback or MHD processes, but underscore 
that self- gravity alone is an important, non-negligible source 
of turbulence in galactic discs. We believe that this work 
is complementary to alternative sources of turb ulence, see 
Sect.[Tl For example. [Hennebelle AuditI (|200/t ) considered 
turbulence driven by colliding fiows in thermally unstable 
gas on very small (parsec) scales which are far from resolved 
in our simulations as we have aimed to resolve the large scale 
contribution from self-gravity that still would be within the 
large beam size {^^ 700 pc). 

Other s tudies of large scale galactic turbulence includes 
IWada et al] (|2002) and Wada & Norman (2007) who used 
an Eulerian code to simulate the dense central part of a 
galactic disc, where the cold molecular gas phase is domi- 
nating. Their results are in agreement with that found here, 
showing a complic ated ISM with a wide range of Q- values. 
Using SPH, Gerri tsen Ickd (I1997I) studied star formation 
and global evolution of the gas in a disc similar to NGC 
650S. They demonstrated that a transient fiocculant spiral 
structure with cold cloud complexes is naturally produced 
in the cold gas, in agreement with our results. The larger 
amount of warm gas was attributed to heating from stel- 
lar phot ons which is ne glected in our work. The subsequent 
work bv lBottemal (|200Sl ) extended parameter space to under- 
stand the relationship between disc mass and global spiral 
structure and pointed out the success of swing amplification 
in predicting this. The measured gas velocity dispersion is 
similar to that obtained in this paper but was attributed to 
mechanical forcing from supernovae feedback. 

Future work attempting more realistic formation of 
molecular clouds requires higher resolution and more so- 
phisticated modelling of radiative physics and feedback in 
order to recover their full range of sizes, masses and life- 
times, which are affected by internal turbulence and strong 
feedback disruption. Even if the actual life and reformation 
times change, we believe that the global evolution of self- 
gravity driven turbulence will remain intact as it is not the 
absolute small scale state of the gas that governs the drag 
of the diffuse gas but the existence of massive interacting 
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agglomerations. These massive clouds form through gravi- 
tational instability that requires seed fluctuations that may 
be triggered initially by numerical noise, but due to their 
rapid growth we expect the long term statistical behaviour 
to be representative. Similarly, our treatment of feedback is 
quite simplistic and could affect the lifetimes of the smaller 
mass clouds at the limit of our resolution, cf. end- state of 
RUNl and RUNS in Fig.[3l although we expect the larger 
clouds to be stable against these effects. Our initial condi- 
tions represent nearby well observed Sc galaxies such as M33 
or NGC 1058. As we initialize the baryonic component as 
gas only we form very massive clouds at early times. How- 
ever, these structure are significantly reduced in mass at late 
times due to star formation. 

It is important to point out that the performed simula- 
tions do not include an old stellar population in the initial 
condition. This might change the global evolution to some 
extent and render a more pronounced spiral structure such 
as an m = 2 mode. Such a setup is rather complex and would 
involve a much larger parameter study. We have postponed 
this to a further study. In addition, we do not include a 
background UV field (far and local field). This will change 
the heating/cooling budget to some extent and can affect 
the gas density distribution. The global dynamics should 
however remain the same (see e.g. IWada et al.ll2002l V 
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